#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mathlink.h"

MLINK link;
MLEnvironment env;

void findfit_(int * num_dim, int * num_monomial, int * num_data, 
              const double * x, const double * y, double * c, const int * exponents)
{
    int packet;

    const int str_size = 1000;
    char model[str_size];
    char param[str_size];
    char var[str_size];
    int i, j, k;

    double * data;
    int dims[2];
    char param_list[str_size];
    const char * text;

    double * res_c;
    int res_n;

    /* ======================================================== */
    /*                     Arguments Part                       */
    /* ======================================================== */

    /* setup parameters */
    k = 0;
    sprintf(&param[k], "{");
    k = k+1;
    for(i = 0; i < *num_monomial-1; ++i) {
        sprintf(&param[k], "c%2.2i,", i+1);
        k = k+4;
    }
    sprintf(&param[k], "c%2.2i}\0", *num_monomial);
    strncpy(param_list, param, strlen(param)+1);
    strcat(param_list, "/.rules");
    /* setup variables */
    k = 0;
    sprintf(&var[k], "{");
    k = k+1;
    for(i = 0; i < *num_dim-1; ++i) {
        sprintf(&var[k], "x%2.2i,", i+1);
        k = k+4;
    }
    sprintf(&var[k], "x%2.2i}\0", *num_dim);
    /* setup fitting model */
    k = 0;
    sprintf(&model[k], "c01");
    k = k+3;
    for(i = 1; i < *num_monomial; ++i) {
        sprintf(&model[k], "+");
        k = k+1;
        sprintf(&model[k], "c%2.2i", i+1);
        k = k+3;
        for(j = 0; j < *num_dim; ++j) {
            sprintf(&model[k], "*x%2.2i^%i", j+1, 
                *(exponents+j*(*num_monomial)+i));
            k = k+6;
        }
    }
    /* setup data */
    data = (double *) calloc((*num_data)*(1+*num_dim), sizeof(double));
    k = 0;
    for(j = 0; j < *num_data; ++j) {
        for(i = 0; i < *num_dim; ++i) {
            *(data+k) = *(x+j*(*num_dim)+i);
            k = k+1;
        }
        *(data+k) = *(y+j);
        k = k+1;
    }
    dims[0] = *num_data;
    dims[1] = *num_dim+1;

    /* ======================================================== */
    /*                      MathLink Part                       */
    /* ======================================================== */

    /* sending FindFit packet to MathLink */
    MLPutFunction(link, "EvaluatePacket", 1);
        MLPutFunction(link, "CompoundExpression", 2);
            MLPutFunction(link, "Set", 2);
                MLPutSymbol(link, "rules");
                MLPutFunction(link, "FindFit", 4);
                    MLPutReal64Array(link, data, dims, NULL, 2);
                    MLPutFunction(link, "ToExpression", 1);
                        MLPutString(link, model);
                    MLPutFunction(link, "ToExpression", 1);
                        MLPutString(link, param);
                    MLPutFunction(link, "ToExpression", 1);
                        MLPutString(link, var);
            MLPutSymbol(link, "Null");
    MLEndPacket(link);

    /* discarding return packets of FindFit */
    while((packet = MLNextPacket(link)) && packet != RETURNPKT)
        MLNewPacket(link);
    MLNewPacket(link);

    /* using return rules of FindFit to get real number list */
    MLPutFunction(link, "EvaluatePacket", 1);
        MLPutFunction(link, "ToExpression", 1);
            MLPutString(link, param_list);
    MLEndPacket(link);

    /* receiving packets from MathLink */
    while((packet = MLNextPacket(link)) && packet != RETURNPKT) {
        printf("packet %i\n", packet);
        if(packet == TEXTPKT) {
            MLGetString(link, &text);
            printf("TEXTPKT:\n");
            printf("%s\n", text);
            MLReleaseString(link, text);
        }
        MLNewPacket(link);
    }
    if(!packet) {
        printf("packets receiving error: %s\n", MLErrorMessage(link));
    } else {
        if(!MLGetReal64List(link, &res_c, &res_n)) {
            printf("MLGetReal64List: internal errror: %s\n", MLErrorMessage(link));
            return;
        }
    }

    for(i = 0; i < res_n; ++i)
        c[i] = res_c[i];

    MLReleaseReal64List(link, res_c, res_n);

    free(data);

    return;
}
